Maestría en Explotación de datos y Descubrimiento de conocimiento


Materia: Enfoque Estadistico del Aprendizaje
Trabajo práctico 1: Regresion Lineal

2. Modelo inicial

Se plantea la siguiente primera alternativa para modelar el peso:

\(E(peso) = \beta_0 + \beta_1 * altura + \beta_2 * edad + \beta_3 * genero + \beta4 * diasActividadFisicaSemanal + \beta5 * consumoDiarioAlcohol\)

Primero se cargan las librerías necesarias:

options(warn=-1)
rm(list=ls())
gc()

         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 535583 28.7    1197378   64   665717 35.6
Vcells 846550  6.5    8388608   64  1812824 13.9

options(warn=-2)

# install.packages("pacman") -- Descomentar par instalar pacman
library(pacman)
p_load_gh('adrianmarino/commons')
import('../src/dataset.R')

[1] "-> '../src/dataset.R' script loadded successfuly!"

import('../src/preprocessing.R')

[1] "-> '../src/preprocessing.R' script loadded successfuly!"

import('../src/model.R')

[1] "-> '../src/model.R' script loadded successfuly!"

import('../src/plot.R')

[1] "-> '../src/plot.R' script loadded successfuly!"

A continuación se carga los conjuntos de entrenamiento y test. también se resumen los valores de las variables categóricas y se excluyen las observaciones con valores faltantes, ya que son muy pocas con redspecto al total.

train_set <- load_train_set() %>% preprocess() %>% shorten_values() %>% drop_missings()

Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(mising_column_names)` instead of `mising_column_names` to silence this message.
ℹ See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
This message is displayed once per session.

test_set  <- load_test_set() %>% preprocess() %>% shorten_values() %>% drop_missings()

glimpse(train_set)

Rows: 6,755
Columns: 15
$ edad                          <int> 17, 15, 15, 16, 17, 15, 13, 17, 17, 16, 16, 14, 15, 17, 15, 14, 17, 17, 16, 14, 1…
$ genero                        <fct> Femenino, Masculino, Masculino, Masculino, Masculino, Masculino, Femenino, Femeni…
$ nivel_educativo               <ord> 2, 1, 2, 1, 2, 1, 9, 9, 1, 3, 3, 8, 9, 3, 9, 2, 3, 3, 2, 9, 8, 2, 3, 2, 2, 3, 1, …
$ altura                        <int> 165, 178, 172, 170, 170, 178, 156, 163, 164, 167, 185, 146, 180, 175, 183, 165, 1…
$ peso                          <int> 62, 62, 62, 65, 75, 88, 46, 60, 57, 51, 100, 33, 62, 70, 80, 60, 50, 50, 70, 75, …
$ frecuencia_hambre_mensual     <ord> Rara vez, Rara vez, Nunca, Nunca, Rara vez, Nunca, Nunca, Nunca, Nunca, Nunca, Nu…
$ dias_consumo_comida_rapida    <int> 0, 0, 3, 1, 1, 2, 0, 0, 0, 3, 4, 2, 1, 1, 3, 0, 0, 0, 1, 0, 6, 0, 1, 0, 2, 0, 2, …
$ edad_consumo_alcohol          <ord> 14-15, <=7, 0, 14-15, 16-17, 8-9, 10-11, 16-17, <=7, 0, 12-13, 12-13, 0, 14-15, <…
$ consumo_diario_alcohol        <dbl> 5.0, 4.0, 0.0, 0.0, 0.0, 5.0, 1.0, 0.5, 5.0, 0.0, 5.0, 0.0, 0.0, 2.0, 1.0, 0.0, 0…
$ dias_actividad_fisica_semanal <int> 7, 7, 7, 7, 0, 7, 0, 2, 7, 3, 2, 2, 7, 1, 4, 0, 6, 5, 7, 3, 0, 7, 5, 2, 2, 4, 2, …
$ consumo_semanal_frutas        <ord> 0, 0, 0, 4-6, 14, 7, 14, 21, 0, 14, <=3, <=3, 7, <=3, <=3, <=3, <=3, <=3, 14, <=3…
$ consumo_semanal_verdura       <ord> 4-6, 4-6, 7, >=28, <=3, 14, 4-6, 7, 0, 4-6, <=3, 7, 7, <=3, 4-6, <=3, <=3, 4-6, <…
$ consumo_semanal_gaseosas      <ord> <=3, <=3, 4-6, <=3, 7, 4-6, 0, 7, <=3, 4-6, 4-6, <=3, <=3, 4-6, 4-6, <=3, 0, 0, 7…
$ consumo_semanal_snacks        <ord> <=3, 0, 4-6, <=3, 0, 4-6, 0, <=3, 0, <=3, <=3, 0, <=3, 7, <=3, 0, <=3, <=3, <=3, …
$ consumo_semanal_comida_grasa  <ord> 0, 4-6, 0, 0, <=3, 4-6, <=3, 7, 0, <=3, 0, <=3, 0, 7, 0, 4-6, 0, <=3, <=3, <=3, <…

Se fija la semilla y se validan las proporciones de los conjuntos de entrenamiento y test:

set.seed(25)
show_train_test_props(train_set, test_set)

[1] "Train: 70%, Test: 30%"

Modelo 1

Se plantea el primer modelo lineal:

model_1 <- lm(
  peso ~ altura + edad + genero + dias_actividad_fisica_semanal + consumo_diario_alcohol, 
  data = train_set
)

¿Cuál es la interpretación de cada uno de los coeficientes estimados?

Veamos a continuación un resumen de los coeficiente del modelo 1:

coefficients_summary(model_1)

__________________________________________________________________________________________________________ 
term                               estimate  std.error   statistic       p.value    conf.low     conf.high 
========================================================================================================== 
(Intercept)                   -68.924840224 2.38253736 -28.9291750 1.413996e-173 -73.5953653 -6.425432e+01 
altura                          0.653666837 0.01466302  44.5792639  0.000000e+00   0.6249227  6.824110e-01 
edad                            1.378761653 0.09535841  14.4587312  1.095370e-46   1.1918291  1.565694e+00 
generoMasculino                 1.224389872 0.27768320   4.4093049  1.053110e-05   0.6800432  1.768737e+00 
dias_actividad_fisica_semanal  -0.099158550 0.05092770  -1.9470456  5.157069e-02  -0.1989929  6.758077e-04 
consumo_diario_alcohol          0.008569358 0.06257711   0.1369408  8.910817e-01  -0.1141015  1.312402e-01 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
____________________________________________________________________________ 
Termino                         Coeficiente Signiticativo IC incluye al cero 
============================================================================ 
(Intercept)                   -68.924840224            Si                 No 
altura                          0.653666837            Si                 No 
edad                            1.378761653            Si                 No 
generoMasculino                 1.224389872            Si                 No 
dias_actividad_fisica_semanal  -0.099158550            No                 Si 
consumo_diario_alcohol          0.008569358            No                 Si 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

Al analizar cada coeficiente se encuentra que:

• \(\hat{\beta_0}\) (Ordenada al origen) de valor -68.92 Kg, es el peso esperado o promedio de un individuo de genero femenino que tiene cero altura, edad, actividad física y consumo diario de alcohol. Esto no es interpretable, ya que una persona tiene que tener una altura superior a cero y no puede tener un peso negativo, pero si podría no realizar actividad física ni consumir alcohol.

• El coeficiente \(\hat{\beta_1}\) de valor 653 gramos, corresponde a la altura del individuo. Este coeficiente indica que dada una edad, genero, consumo de alcohol diario y días de actividad física semanal fijos, cada incremento en 1 cm adicional en la altura del individuo implica un aumento de su peso esperado o promedio de 653 gramos.

• El coeficiente \(\hat{\beta_2}\) de valor 1.378 kg, corresponde a la edad del individuo. Este coeficiente indica que dada una altura, genero, días de actividad física y consumo de alcohol diario fijos, cada vez que el individuo cumple un año su peso esperado o promedio aumenta en 1.378 kg.

• El coeficiente \(\hat{\beta_3}\) de valor 1.224 kg, corresponde a los individuos de genero masculinos. Este coeficiente indica que dada una altura, edad, consumo de alcohol diario y días de actividad física semanal fijos, el peso promedio o esperado para el genero masculino es 1.224 kg mayor al peso femenino (categoría basal). Por otro lado, el coeficientes no indica cunado mas alto es el peso del genero masculino respecto del femenino al fijar los demás coeficientes.

• El coeficiente \(\hat{\beta_4}\) de valor 99.1 gramos, corresponde a los días de actividad física semanal que realiza el individuo. Este coeficiente indica que dada una altura, edad, genero y consumo de alcohol diario, cada vez que un individuo realiza un día mas de actividad física semanal su peso esperado o promedio disminuye en 99.1 gramos.

• El coeficiente \(\hat{\beta_5}\) de valor -8 gramos, corresponde al nivel de consumo diario de alcohol del individuo. Este coeficiente indica que dada una altura, edad, genero y días de actividad física semanal fijos, cada vez que el individuo consume un trago de alcohol su peso esperado o promedio disminuye en 8 gramos. A simple vista podrá no llegar a tener sentido, ya que a mayor consumo de alcohol el peso debería aumentar, ya sea por el peso del propio liquido como el peso equivalente en grasas. Entiendo que puede tener un relación con los rangos de edades de los individuos que mas consumen alcohol (12 q 17 años), ya que estos se encuentran en pleno crecimiento.

¿Son significativos los coeficientes?

Para determina si los coeficientes son aptos para explicar el peso de un individuo se realiza un \({T}\) test para cada coeficiente en el cual se evalúan las siguientes hipótesis:

• \({H_0: \beta_i = 0}\)
• \({H_1: \beta_i \neq 0}\)

Si \({\beta_i \neq 0}\) podemos decir que existe una diferencia estadisticamente significativas del cero para coeficiente \({\beta_i}\), y por lo tanto el coeficiente \({\beta_i}\) explicar la variable \({y}\) (Peso en nuestro caso).

Luego analizando la salida de coefficients_summary concluimos que:

• Los coeficientes correspondientes al acategoria basal \({\beta_1}\)(Genero femenino), altura(\({\beta_1}\)), edad(\({\beta_2}\)) y genero masculino (\({\beta_3}\)) tienen \(p-valor < 0.05\). Por lo tanto, se rechaza la hipótesis nula y resultan estadistitamente significativos para explicar el peso. Por otro lado, se puede apreciar que los intervalos de confianza del 95% no incluyen al cero.
• Lo contrario sucede con días de actividad física semanal(\({\beta_4}\)) y consumo de alcohol diario (\({\beta_5}\)), dado que ambos no rechazar la hipótesis nula (\(p-valor > 0.05\)) y por lo tanto no existe una diferencia significativa del cero. Finalmente, no hay evidencia estadistitamente significativas de que estos coeficientes expliquen al peso.

¿El modelo resulta significativo para explicar el peso?

Para determinar si es modelo es significativo para explicar el peso de un individuo se realiza un \(F\) test con las siguientes hipótesis:

• \(H_0: β_1 = β_2 = · · · = β_{p−1} = 0\)
• \(H_1:\) Por lo menos un \(β_k\) (\(k = 1, 2,..., p−1\)) es distinto de 0.

Donde: * \(H_0\) afirma que no hay vinculo entre la variable \({y}\)(Peso) y las variables regresoras. * \(H_1\) afirma que al menos una de las variables regresoras sirve para predecir la variable \({y}\) (Peso).

Veamos los resultados el \(F\) test:

glance(model_1)

Podemos apreciar que el \(p-valor < 0.05\) e igual a 0. Con mucha certeza podemos decir que al menos una de las variables regresoras permite explicar el peso. Esto concuerda con los resultados de los \(T\) test para las los coeficientes correspondientes a altura, edad y genero femenino(basa) y masculino).

¿Qué porcentaje de la variabilidad explica el modelo?

Según el valor de \(R^2\) ajustado (adj.r.squared), este modelo llega a explica el 35% de la variabilidad del dataset de entrenamiento, lo cual no es un valor bajo pero tampoco es despreciable.

¿Que sucede si poner al genero masculino como variable basal?

train_set_genero <- data.frame(train_set) 
train_set_genero$genero <- factor(
  train_set_genero$genero,
  levels=c('Masculino', 'Femenino'), 
  ordered=FALSE
)
table(train_set_genero$genero)

Masculino  Femenino 
     3100      3655 

model_genero <- lm(
  peso ~ altura + edad + genero + dias_actividad_fisica_semanal + consumo_diario_alcohol, 
  data = train_set_genero
)

coefficients_summary(model_genero)

__________________________________________________________________________________________________________ 
term                               estimate  std.error   statistic       p.value    conf.low     conf.high 
========================================================================================================== 
(Intercept)                   -67.700450352 2.49636490 -27.1196131 7.976726e-154 -72.5941133 -6.280679e+01 
altura                          0.653666837 0.01466302  44.5792639  0.000000e+00   0.6249227  6.824110e-01 
edad                            1.378761653 0.09535841  14.4587312  1.095370e-46   1.1918291  1.565694e+00 
generoFemenino                 -1.224389872 0.27768320  -4.4093049  1.053110e-05  -1.7687366 -6.800432e-01 
dias_actividad_fisica_semanal  -0.099158550 0.05092770  -1.9470456  5.157069e-02  -0.1989929  6.758077e-04 
consumo_diario_alcohol          0.008569358 0.06257711   0.1369408  8.910817e-01  -0.1141015  1.312402e-01 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
____________________________________________________________________________ 
Termino                         Coeficiente Signiticativo IC incluye al cero 
============================================================================ 
(Intercept)                   -67.700450352            Si                 No 
altura                          0.653666837            Si                 No 
edad                            1.378761653            Si                 No 
generoFemenino                 -1.224389872            Si                 No 
dias_actividad_fisica_semanal  -0.099158550            No                 Si 
consumo_diario_alcohol          0.008569358            No                 Si 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

glance(model_genero)

Observaciones

• En este caso, cambia el coeficiente \(\β_0\) para la categoría basal, manteniéndose los demás coeficientes.
• Siguen siendo significativos los mismo coeficientes. La unica diferencia es que \(\β_0\)(basal) representa al genero masculino y \(\β_3\) al femenino.
• Se mantiene la significatividad global del modelo y el mismo \(R^2\).
• Como ultimo, cabe aclarar que al cambiar la categoría basar cambian las interpretaciones del modelo.

3. Modelo categóricas

Se sugiere probar un modelo que incorpore el consumo semanal de snacks y una interacción entre el género y la edad, en lugar de actividad física y consumo de alcohol. Además se pide explicitamente que la categoría “No comí comida salada o snacks en los últimos 7 días” de la variable consumo_semanal_snacks se encuentre como nivel/categoría basal.

\(E(peso) = \beta_0 + \beta_1 * altura + \beta_2 * edad + \beta_3 * genero + \beta4 * consumoSemanalSnacks + \beta_5 * genero * edad\)

Primero validamos que las primeras categorías en cada variable de tipo factor sean las correctas, ya que esta sera la que el modelo defina como categoría basal:

table(train_set$consumo_semanal_snacks)

   0  <=3  4-6    7   14   21 >=28 
2090 3044  599  580  223   96  123 

table(train_set$genero)

 Femenino Masculino 
     3655      3100 

Se puede apreciar que la primeras categorías corresponden a 0 consumo de snacks semanal y genero femenino. Por otro lado la categoría genero se encuentra balanceada.

Modelo 2

Definimos el nuevo modelo:

model_2 <- lm(
  peso ~ altura + edad + genero + consumo_semanal_snacks +  genero * edad, 
  data = train_set
)

¿Cuál es la interpretación de los coeficientes estimados para las categorías de consumo_semanal_snacks y genero*edad? ¿Son significativas?

coefficients_summary(model_2)

____________________________________________________________________________________________________ 
term                          estimate  std.error  statistic       p.value      conf.low   conf.high 
==================================================================================================== 
(Intercept)                -64.6479875 2.87030372 -22.523048 2.261394e-108 -70.274689295 -59.0212858 
altura                       0.6467154 0.01481585  43.650234  0.000000e+00   0.617671609   0.6757591 
edad                         1.2167591 0.12239655   9.941123  3.974656e-23   0.976823190   1.4566950 
generoMasculino             -4.0403437 2.72370910  -1.483398  1.380155e-01  -9.379673695   1.2989863 
consumo_semanal_snacks<=3   -1.4368851 0.28028645  -5.126488  3.034134e-07  -1.986335069  -0.8874352 
consumo_semanal_snacks4-6   -2.2551628 0.45726634  -4.931837  8.340962e-07  -3.151549278  -1.3587764 
consumo_semanal_snacks7     -0.6565845 0.46383525  -1.415555  1.569519e-01  -1.565848025   0.2526791 
consumo_semanal_snacks14    -1.0608560 0.69485004  -1.526741  1.268724e-01  -2.422981482   0.3012696 
consumo_semanal_snacks21    -1.1903835 1.02955288  -1.156214  2.476346e-01  -3.208632258   0.8278653 
consumo_semanal_snacks>=28  -2.6018645 0.91604489  -2.840324  4.520270e-03  -4.397601785  -0.8061272 
edad:generoMasculino         0.3493851 0.18203251   1.919355  5.498154e-02  -0.007456148   0.7062263 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
_______________________________________________________________________ 
Termino                    Coeficiente Signiticativo IC incluye al cero 
======================================================================= 
(Intercept)                -64.6479875            Si                 No 
altura                       0.6467154            Si                 No 
edad                         1.2167591            Si                 No 
generoMasculino             -4.0403437            No                 Si 
consumo_semanal_snacks<=3   -1.4368851            Si                 No 
consumo_semanal_snacks4-6   -2.2551628            Si                 No 
consumo_semanal_snacks7     -0.6565845            No                 Si 
consumo_semanal_snacks14    -1.0608560            No                 Si 
consumo_semanal_snacks21    -1.1903835            No                 Si 
consumo_semanal_snacks>=28  -2.6018645            Si                 No 
edad:generoMasculino         0.3493851            No                 Si 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

Si interpretamos los coeficientes que son significativos para el \(T\) test:

• Coeficiente correspondiente al consumo_semanal_snacks<=3: Si fijamos los coeficientes correspondientes a la altura, edad, generoMasculino y generoMasculino*edad; el peso promedio o esperado de un individuo de consume snacks hasta 3 veces por semana es 1.43 kg menor que aquellos que no consumen snacks.

• Sucede algo similar con las categorias consumo_semanal_snacks4-6 y consumo_semanal_snacks>=28 donde:

  • El peso para un individio que consume de 4 a 6 veces por semana es 2.25 kg menor que aquellos con no consumen snacks.
  • El peso para un individio que consume 28 veces a mas por semana es 2.6 kg menor que aquellos con no consumen snacks.

• Coeficiente correspondiente al edad:generoMasculino: Si fijamos los coeficientes correspondiente a la altura, edad, genero masculino y consumo de snack; el peso promedio o esperado de un individuo de genero masculino por edad es 0.349 kg mayor a la categoria basal edad:generoFemenino.

¿Qué porcentaje de la variabilidad explica el modelo? En caso de detectar que existen categorías no significativas de la variable consumo_semanal_snacks evaluar si la variable es significativa en su conjunto y, en caso afirmativo, proponer una redefinición de las mismas que permita obtener una mayor proporción de categorías significativas individualmente. Luego, analizar si existen cambios en la variabilidad explicada por el modelo.

Viendo el resultado de coefficients_summary se aprecia que las siguientes categorías de consumo_semanal_snacks no son significativas:

• 2 veces al día (14 veces/semana)
• 4 a 6 veces durante los últimos 7 (4 a 6 veces semana)
• 3 veces al día (21 veces/semana)

Pero si son significativas los extremos:

• De 1 a 3 veces/semana
• De 28 o mas veces/semana

A continuación se realiza un \(F\) test para evaluar la significatividad conjunta de las categóricas de la variable consumo_semanal_snacks para explicar el peso.

El \(F\) test también llamando ANOVA (Análisis de la varianza) se realiza para probar la significatividad conjunta de todos los valores de una variable categórica.

Las hipótesis son las siguientes:

• \(H_0: β_q = β_{q+1} = · · · = β_{p−1} = 0\)

• \(H_1:\) por lo menos uno de los \(β_k\) (con \(k\) entre \(q\) y \(p−1\)) es tal que \(β_k \neq 0\).

Luego si todos los coeficientes asociados a los valores de variable categórica son cero, se rechaza la hipótesis nula y por lo tanto la variable no es significartiva para explicar el peso en nuestro caso.

A continuación veremos el p-valor resultado de aplicar \(F\) test para cada variable del modelo:

anova_summary(model_2)

Podemos apreciar que el \(p-value < 0.005\) para la variable consumo_semanal_snacks. Por lo tanto se rechaza la hipótesis nula y podemos decir en su conjunto resulta estadísticamente significativa para explicar el peso. Luego, como la variable consumo_semanal_snacks es significativa vale la pena re-definirla. Por otro lado, la combinación de variables genero-edad no es estadísticamente significativa para explicar el peso, pero si lo es el genero en forma separada. Finalmente, como ya vimos en pasos anteriores, edad y altura son significativas.

Veamos a continuación las distribuciones de las categorías de la variable consumo_semanal_snacks ordenadas por la mediana del peso:

segmented_box_plot(
  train_set, 
  column        = 'peso', 
  segmented_by  = 'consumo_semanal_snacks',
  title         = 'Consumo de snacks ordenado por la mediana del peso',
  y_label       = 'Peso (Kg)',
  y_limits      = c(10, 130),
  x_label       = 'Consumo de snacks (Veces/Semana)'
)

A simple vista no parece haber una gran diferencia, pero si se aprecia que los extremos difieren del los valores centrales.

A continua con se promane una nueva definición de la variable consumo_semanal_snacks. Primero se realiza el promediodel peso para cada categoría de la variable consumo_semanal_snacks:

peso_medio_by_nivel_consumo_snack = train_set %>% 
  group_by(consumo_semanal_snacks) %>%
  summarise(promedio = mean(peso))

ggplot(data = peso_medio_by_nivel_consumo_snack, aes(x = promedio)) + 
  geom_boxplot(alpha = 0.75, fill="blue") +
  labs(title = "Peso promedio por cada categoria de consumo de snacks") +
  labs(x = "Peso medio") +
  theme_bw()

peso_medio_by_nivel_consumo_snack  <- peso_medio_by_nivel_consumo_snack$promedio

peso_medio_by_nivel_consumo_snack

[1] 60.68182 59.13633 58.29215 58.29483 58.69507 58.27083 55.52033

quantile(peso_medio_by_nivel_consumo_snack)

      0%      25%      50%      75%     100% 
55.52033 58.28149 58.29483 58.91570 60.68182 

Se puede apreciar que es una distribución asimétrica sesgada a derecha, ya que los mayores valores se encueran arriba del segundo cuantil (Mediana).

A continuación se re-definen las categorías originales por 3 nueva categorías: Bajo, Medio, Alto. Esta categorías estan asociadas al peso de de individuo. Si el indivídio tiene un peso menor al Q1, se le asigna el nivel Bajo, si esta entre Q1 y Q3 sera Medio y Alta arroba de Q3. Finalmente a continuación se transforma la variable en el conjunto de entrenamiento y en test se usas los los cuantíles generados con en conjunto de entrenamiento:

q1 <- quantile(peso_medio_by_nivel_consumo_snack)[2]
q3 <- quantile(peso_medio_by_nivel_consumo_snack)[4]

train_set2 <- train_set %>%
  mutate(consumo_semanal_snacks = case_when(peso <  q1 ~ "Bajo", peso >= q3 ~ "Alto", TRUE ~ "Medio"))

test_set2 <- test_set %>%
  mutate(consumo_semanal_snacks = case_when(peso <  q1 ~ "Bajo", peso >= q3 ~ "Alto", TRUE ~ "Medio")) %>% 
  mutate(consumo_semanal_snacks = as.factor(consumo_semanal_snacks))

test_set2 %>% segmented_box_plot(
  column        = 'peso', 
  segmented_by  = 'consumo_semanal_snacks',
  title         = 'Consumo de snacks ordenado por la mediana del peso en Test',
  y_label       = 'Peso (Kg)',
  y_limits      = c(40, 100),
  x_label       = 'Consumo de snacks (Veces/Semana)'
)

Se puede apreciar que en el conjunto de test no hay valores medio, pero existe en el conjunto de entrenamiento.

Modelo 3

A continuación definimos un nuevo modelo igual al anterior pero ahora ya usando la re-definición de la variable consumo_semanal_snacks:

model_3 <- lm(
  peso ~ altura + edad + genero + consumo_semanal_snacks +  genero * edad, 
  data = train_set2
)

models <- list('Modelo 1'=model_1, 'Modelo 2'=model_2, 'Modelo 3'=model_3)

coefficients_summary(model_3)

__________________________________________________________________________________________________ 
term                          estimate  std.error  statistic       p.value    conf.low   conf.high 
================================================================================================== 
(Intercept)                  9.1095905 2.33677462   3.898361  9.778446e-05   4.5287749  13.6904061 
altura                       0.3026553 0.01186350  25.511464 3.919595e-137   0.2793991   0.3259115 
edad                         0.5805838 0.09095052   6.383513  1.845540e-10   0.4022921   0.7588755 
generoMasculino             -2.2818007 2.01444847  -1.132717  2.573732e-01  -6.2307553   1.6671540 
consumo_semanal_snacksBajo -15.5664145 0.20735911 -75.069836  0.000000e+00 -15.9729038 -15.1599252 
edad:generoMasculino         0.1631882 0.13464886   1.211954  2.255725e-01  -0.1007661   0.4271424 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
_______________________________________________________________________ 
Termino                    Coeficiente Signiticativo IC incluye al cero 
======================================================================= 
(Intercept)                  9.1095905            Si                 No 
altura                       0.3026553            Si                 No 
edad                         0.5805838            Si                 No 
generoMasculino             -2.2818007            No                 Si 
consumo_semanal_snacksBajo -15.5664145            Si                 No 
edad:generoMasculino         0.1631882            No                 Si 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

anova_summary(model_3)

models %>% 
  map_df(glance, .id = "model") %>%
  arrange(desc(adj.r.squared))

Finalmente se aprecia que la nueva categorización de la variable aumenta el \(R^2\) Ajustado a casi el doble (Del 35.7 al 64.7%). Si bien esto mejora la capacidad explicativa del modelo, en pasos posteriores se deberá determina si produce o no overfitting sobre el conjunto de entrenamiento.

4. Modelos propios y evaluación

Realizar 2 modelos lineales múltiples adicionales y explicar breve-mente la lógica detrás de los mismos (se valorará la creación y/o inclusión de variables nuevas).

Evaluar la performance del modelo inicial, el modelo categóricas con las categorías redefinidas de la variable consumo_semanal_snacks y los modelos desarrollados en este punto en el dataset de entrenamiento y evaluación (usar dataset “encuesta_salud_test.csv”).

La evaluación de performance consiste en comparar en ambos sets la performance en términos del R cuadrado ajustado, RMSE y MAE.

Al continuación se define 2 modelos.

Modelo 4

\(E(peso) = \beta_0 + \beta_1 * altura + \beta_2 * edad + + \beta_3 * genero + \beta4 * consumoSemanalSnacks + \beta_5 * diasActividadFisicaSemanal + \beta_6 * altura * genero\)

Se utilizo la redefinición de la variable consumo_semanal_snacks como base. Ademase se agregar la variable dias_actividad_fisica_semanal entendiendo que tiene una influencia iportante en el peso y luego la asociacion altura * genero ya que en general mas mujeres tienen a ser mas bajar que los varones y vise versa.

model_4 <- lm(
  peso~ 
    altura + 
    edad + 
    genero + 
    consumo_semanal_snacks +
    dias_actividad_fisica_semanal +
    altura*genero,
  data = train_set2
)

coefficients_summary(model_4)

______________________________________________________________________________________________________ 
term                             estimate  std.error  statistic      p.value     conf.low    conf.high 
====================================================================================================== 
(Intercept)                    21.3416310 2.92913011   7.285996 3.556273e-13  15.59961152  27.08365041 
altura                          0.2259395 0.01724260  13.103560 9.312479e-39   0.19213856   0.25974045 
edad                            0.6053145 0.06966453   8.688990 4.534512e-18   0.46874999   0.74187894 
generoMasculino               -22.1352833 3.47288831  -6.373739 1.966410e-10 -28.94324039 -15.32732614 
consumo_semanal_snacksBajo    -15.5990623 0.20665882 -75.482200 0.000000e+00 -16.00417882 -15.19394580 
dias_actividad_fisica_semanal  -0.1055494 0.03741373  -2.821141 4.799296e-03  -0.17889212  -0.03220668 
altura:generoMasculino          0.1364081 0.02117342   6.442424 1.256771e-10   0.09490156   0.17791473 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
__________________________________________________________________________ 
Termino                       Coeficiente Signiticativo IC incluye al cero 
========================================================================== 
(Intercept)                    21.3416310            Si                 No 
altura                          0.2259395            Si                 No 
edad                            0.6053145            Si                 No 
generoMasculino               -22.1352833            Si                 No 
consumo_semanal_snacksBajo    -15.5990623            Si                 No 
dias_actividad_fisica_semanal  -0.1055494            Si                 No 
altura:generoMasculino          0.1364081            Si                 No 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

anova_summary(model_4)

glance(model_4)

train_set3 <- column_mean_quantile_binning(train_set2, 'dias_actividad_fisica_semanal')
test_set3  <- column_mean_quantile_binning(test_set2,  'dias_actividad_fisica_semanal')

train_set3 <- column_mean_quantile_binning(train_set3, 'consumo_semanal_frutas')
test_set3  <- column_mean_quantile_binning(test_set3,  'consumo_semanal_frutas')

train_set3 <- column_mean_quantile_binning(train_set3, 'consumo_semanal_verdura')
test_set3  <- column_mean_quantile_binning(test_set3,  'consumo_semanal_verdura')

train_set3 <- column_mean_quantile_binning(train_set3, 'consumo_semanal_comida_grasa')
test_set3  <- column_mean_quantile_binning(test_set3,  'consumo_semanal_comida_grasa')

train_set3 <- column_mean_quantile_binning(train_set3, 'consumo_semanal_gaseosas')
test_set3  <- column_mean_quantile_binning(test_set3,  'consumo_semanal_gaseosas')

segmented_box_plot(
  test_set3, 
  column        = 'peso', 
  segmented_by  = 'dias_actividad_fisica_semanal',
  title         = 'Niveles actividad fisica ordenados por la mediana del peso en Test',
  y_label       = 'Peso (Kg)',
  y_limits      = c(40, 100),
  x_label       = 'Niveles de actividad física (Dias)'
)

Modelo 5

\(E(peso) = \beta_0 + \beta_1 * altura + \beta_2 * edad + + \beta_3 * genero + \beta4 * consumoSemanalSnacks + \\ \beta_5 * diasActividadFisicaSemanal + \beta_6 * consumoSemanalFrutas + \beta_7 * consumoSemanalVerduras + \\* \beta_8 * consumoSemanalGrasas + \beta_9 * consumoSemanalGaseosas\)

Se utilizo la redefinición de la variable consumo_semanal_snacks como base. Ademase se agregar la variable consumo_semenal_frutras/verduras/grasas/gaseaosas entendiendo que también tiene una influencia importante en el peso.

model_5 <- lm(
  peso ~ 
  edad +
  genero +
  altura +
  consumo_semanal_snacks +
  consumo_semanal_frutas + 
  consumo_semanal_verdura,
  data = train_set3
)

coefficients_summary(model_5)

______________________________________________________________________________________________________ 
term                             estimate  std.error   statistic       p.value    conf.low   conf.high 
====================================================================================================== 
(Intercept)                   10.44337415 2.00375738   5.2118955  1.924670e-07   6.5153774  14.3713709 
edad                           0.60102954 0.06872887   8.7449355  2.782235e-18   0.4662993   0.7357598 
generoMasculino                0.05600315 0.20151595   0.2779093  7.810905e-01  -0.3390317   0.4510380 
altura                         0.29569466 0.01157806  25.5392162 2.052362e-137   0.2729980   0.3183913 
consumo_semanal_snacksBajo   -16.08262471 0.20786131 -77.3719010  0.000000e+00 -16.4900985 -15.6751510 
consumo_semanal_frutasBajo             NA         NA          NA            NA          NA          NA 
consumo_semanal_frutasMedio   -8.88749099 0.64594156 -13.7589707  1.679091e-42 -10.1537403  -7.6212417 
consumo_semanal_verduraBajo            NA         NA          NA            NA          NA          NA 
consumo_semanal_verduraMedio           NA         NA          NA            NA          NA          NA 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
__________________________________________________________________________ 
Termino                       Coeficiente Signiticativo IC incluye al cero 
========================================================================== 
(Intercept)                   10.44337415            Si                 No 
edad                           0.60102954            Si                 No 
generoMasculino                0.05600315            No                 Si 
altura                         0.29569466            Si                 No 
consumo_semanal_snacksBajo   -16.08262471            Si                 No 
consumo_semanal_frutasBajo             NA            NA                 NA 
consumo_semanal_frutasMedio   -8.88749099            Si                 No 
consumo_semanal_verduraBajo            NA            NA                 NA 
consumo_semanal_verduraMedio           NA            NA                 NA 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

anova_summary(model_5)

glance(model_5)

c(models, list('Modelo 4'=model_4, 'Modelo 5'=model_5)) %>% 
  map_df(glance, .id = "model") %>%
  arrange(desc(adj.r.squared))

Finalmente, si comparamos los modelos por \(R^2\) Ajustado, se puede apreciar que el modelo 5 (con todas las variables categóricas re-definidas) llega a captar la mayor varianza explicada sobre el dataset de entrenamiento. Por supuesto esto no dice nada acerca de la performance del modelo en test, pero si que tiene la mejor capacidad para extraer información de los dato de entrenamiento.

¿Cuál es el mejor modelo para nuestro objetivo de predecir el peso? ¿Por qué?

Ahora comparamos la performance de todo los modelos al evaluar el error delos mismo al predecir el peso en el conjunto de train y test tanto para RMSE como MAE:

RMSE

custom_models_evaluation_summary(
  model_1, model_2, model_3, model_4, model_5,
  test_set, test_set2, test_set3,
  metric_fn = rmse
)

Si utilizamos la métrica RMSE podemos ver que el modelo 5 tiene el menor error en el conjunto de test. Por otro lados el que tiene la mayor diferencia de error entre test y entrenamiento. Esto nos dice que podría estar sobre-ajustandose al conjunto de entrenamiento. El modelo 3 tiene un error en test muy cercano y ademas tiene un diferencia entre test y train mucho menor. por esto ultimo parece ser el mejor modelo ya que tiene prácticamente el menor error posible y también el menor sobre-ajuste al conjunto de entrenamiento.

MAE

custom_models_evaluation_summary(
  model_1, model_2, model_3, model_4, model_5,
  test_set, test_set2, test_set3,
  metric_fn = mae
)

Si medimos a partir del MAE sucede algo muy similar, El modelo 3 es es que tiene menor error y ademas menos sobre-ajuste.

Finalmente, según ambas metricas el moejor modelo es el Modelo 3.

5. Diagnóstico del modelo

Analizar en profundidad el cumplimiento de los supuestos del modelo lineal para el modelo inicial.

plot(model_1)

Homosedastisidad

Al visualizar el primer gráfico (Residuos vs. Valores ajustados) se puede apreciar que no hay presencia de homocedastrisisdad, ya que los valores predicho, la variabilidad o amplitud de los residuos parece mantenerse con cierta regularidad. Dadas esta condiciones podemos decir que se cumple el supuesto de varianza constante.

Normalidad

Al visualizar el diagrama QQ-Plot podemos observas que en el extremo derecha, el modelo sobre estima el peso del los individuos ya que hay una gran diferencia positiva entre el valor predicho y el valor esperado teórico. lo mis sucede a izquierca pero en menor medida, donde el modelo subestima el valor de peso en comparación al valor esperado teórico. Como dato adiciona este grafito corresponde a una distribución sesgada a derecha, también conocido como sesgo positivo. Finalmente el QQ-Plot no muestra un grado de alejamiento pronunciado de una districion normal teórica y decimos que no se cumple el supuesto de normalidad del modelo.

Apalancamiento (Leverage)

Si observamos el gráfico de Residuos vs Apalacamiento vemos que varias observaciones o individuos que se alejan del cumulo de principal. Estos ejercen un apuntalamiento sobre el valores predicho del modelo a partir de un apalancamiento(leverage) 0.0020 y es mas pronunciado desde 0.0025. Finalmente vemos un grado importante de desvió de las predicciones vs su vor esperado.

A continuación se pueden ver lo individuos que producen mayor apalancamiento(leverage) y por ende sesgo en al predicción del modelo:

augment(model_1) %>%
  filter(.hat>0.00245) %>%
  arrange(.hat)

6. Modelo Robusto

Leer el archivo “encuesta_salud_modelo6.csv”. Este último consiste en el dataset original de train con la incorporación de algunas observaciones adicionales que pueden incluir valores atípicos. En particular, observar la relación entre peso y altura ¿Qué ocurre con estos nuevos datos? Entrenar el modelo inicial con estos nuevos datos y comentar qué se observa en los coeficientes estimados y las métricas de evaluación (R cuadrado ajustado, RMSE y MAE) respecto al modelo entrenado con el set de entrenamiento original. Entrenar un modelo robusto con la misma especificación que el modelo inicial sobre los nuevos datos. Comparar los coeficientes y su performance (RMSE y MAE) respecto al modelo inicial no robusto entrenado en este punto. ¿Qué puede concluir al respecto?

Se carga el conjunto de entrenamiento en crudo,e s decir sin pre-procesamiento. Luego se resumen los valores de las variables categóricas y se eliminan missing values, ya que siguen siendo muy poco casos:

original_train_set <- shorten_values(preprocess(load_original_train_set()))
missings_summary(original_train_set)

new_train_set <- drop_missings(original_train_set)
missings_summary(new_train_set)

nrow(original_train_set)

[1] 7060

nrow(new_train_set)

[1] 6789

Comparemos las distribuciones del peso vs altura en ambos conjunto de entrenamiento:

box_plots(
  train_set %>% select(peso, altura), 
  title = 'Comparativas de distribuciones del peso y la altura'
)

box_plots(
  new_train_set %>% select(peso, altura), 
  title = 'Comparativas de distribuciones del peso y la altura'
)

En el dataset de entrenamiento original la variable peso tiene prácticamente el doble de outliers que el dataset procesado.

Modelo 6

Definimos un modelo igual al modelo 1 pero entrenando en el dataset de entrenamiento original.

model_6 <- lm(
  peso ~ altura + edad + dias_actividad_fisica_semanal + consumo_diario_alcohol, 
  data = new_train_set
)

coefficients_summary(model_6)

________________________________________________________________________________________________________ 
term                              estimate  std.error   statistic       p.value    conf.low    conf.high 
======================================================================================================== 
(Intercept)                   -73.01459357 2.68343292 -27.2093977 8.053958e-155 -78.2749640 -67.75422317 
altura                          0.69037327 0.01560204  44.2489243  0.000000e+00   0.6597884   0.72095816 
edad                            1.31717336 0.11443750  11.5099803  2.246611e-30   1.0928400   1.54150676 
dias_actividad_fisica_semanal  -0.09224986 0.06116426  -1.5082314  1.315419e-01  -0.2121510   0.02765128 
consumo_diario_alcohol          0.01109380 0.07546845   0.1469992  8.831371e-01  -0.1368480   0.15903564 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
___________________________________________________________________________ 
Termino                        Coeficiente Signiticativo IC incluye al cero 
=========================================================================== 
(Intercept)                   -73.01459357            Si                 No 
altura                          0.69037327            Si                 No 
edad                            1.31717336            Si                 No 
dias_actividad_fisica_semanal  -0.09224986            No                 Si 
consumo_diario_alcohol          0.01109380            No                 Si 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 

anova_summary(model_6)

glance(model_6)

print(paste('Disminicion de adj.r.squared:', abs(0.352113 - 0.2734821) * 100, '%'))

[1] "Disminicion de adj.r.squared: 7.86309 %"

Dada la presencia de outliers en la variable peso, el \(R^2\) Ajustado baja con respecto al modelo 1.

models <- list('Modelo 6'=model_6)

models_evaluation_summary(models, train_set, metric_fn = rmse)

models_evaluation_summary(models, train_set, metric_fn = mae)

Por otro lado, aumento el error de predicción tanto en train como en test. Finalmente, el modelo tiene un grado de overfitting mucho mayor que los modelos anteriores, ya que la métrica de evaluación en test y train tiene una diferencia muy pronunciada de 1.7 puntos.

Modelo 7

Definimos un modelo igual al modelo 1 entrenando en el dataset de entrenamiento original y usamos un modelo lineal robusto.

model_7 <- rlm(
  peso ~ altura + edad + dias_actividad_fisica_semanal + consumo_diario_alcohol, 
  data = new_train_set
)
coefficients_summary(model_7)

___________________________________________________________________________________________ 
term                              estimate  std.error   statistic     conf.low    conf.high 
=========================================================================================== 
(Intercept)                   -69.44476006 2.04036494 -34.0354604 -73.44380185 -65.44571827 
altura                          0.66723355 0.01186310  56.2444292   0.64398229   0.69048481 
edad                            1.23827826 0.08701326  14.2309137   1.06773540   1.40882113 
dias_actividad_fisica_semanal  -0.02163933 0.04650663  -0.4652957  -0.11279064   0.06951198 
consumo_diario_alcohol          0.03595659 0.05738291   0.6266080  -0.07651185   0.14842502 
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
[1] "WARN: p.value column is required to make model coefficients summary!\n"
[1] "WARN: p.value column is required to plot tidy coefficients!\n"
NULL

anova_summary(model_7)

models <- list('Modelo 6'=model_6, 'Modelo 7'=model_7)

models_evaluation_summary(models, test_set, metric_fn = rmse)

models_evaluation_summary(models, test_set, metric_fn = mae)

El modelo lineal robusto (Modelo 7) parece tener un menor error de entrenamiento muy cercano al modelo 6, pero tiene mayor sobre- ajuste que el modelo 6, aunque es una diferencia muy baja.

Dado esto, seria una buena selecciono elegir el modelo 7, ya que el sobre ajuste practicamente no cambia y obtenemos un error de predicción en test ligeramente menor.